/*  WECC Interstate Reallocation with bootstrapped damages.

NOTES: All directories and paths should be set in the "Create_TableA5.do" file.

*/

clear all


import delimited `"${data_path}/Generation_and_damages_withCO2_2019-11-24_eGrid_SR.csv"', case(preserve) 

merge 1:1 zip using `"${data_intermediate}/WECCBootsFinal11232019_wide.dta"'
keep if _merge==3
drop _merge
rename zip zipcode
rename STATE state


merge 1:1 zipcode using `"${data_path}/HousingUnits"'
keep if _merge==3
*need to drop zipcodes with missing housing units or zip codes that have 0 housing units
*Note: There are 121 zip codes with a positive number of installations but 0 housing units
*This can be observed by uncommenting the next line.
*count if single_units ==0 & num_systems>0 & !missing(num_systems)
drop if missing(single_units) | single_units==0
drop _merge

foreach var of varlist total_damages1b1-total_damages8b99 {
	gen total_damages_all = `var'*num_systems
	egen SumOld`var'=sum(total_damages_all)
	drop total_damages_all

}
gen total_damages_all = total_damages*num_systems
egen sumOldPtEst = sum(total_damages_all)
drop total_damages_all

egen old_p25 = rowpctile( SumOldtotal_damages1b1- SumOldtotal_damages8b99), p(2.5)
egen old_p975 = rowpctile( SumOldtotal_damages1b1- SumOldtotal_damages8b99), p(97.5)
egen old_mean = rowmean( SumOldtotal_damages1b1- SumOldtotal_damages8b99)

***PAUSE HERE to record estimates and CIs of avoided damages from currently installed capacity

*drop SumOldtotal_damages1b1- SumOldtotal_damages8b99
*drop old_p25 old_p975 old_mean sumOldPtEst

***Get new damages based on point estimate
gen total_damages_all = total_damages*num_systems
	sort total_damages single_units
	egen rank=rank(-total_damages)
	sort rank single_units
	drop rank
	gen rank=_n /*ranks zipcodes by total_damages and number of houses*/


* Simulation 1 (moving all panels to highest env value zipcode subject to #of housing units)
gen new_num_systems=0

egen sumsystems = sum(num_systems)
gen remaining = sumsystems
sum rank, meanonly
local maxindex=r(max)
forvalues i=1/`maxindex' {
	quietly replace new_num_systems=min(remaining[`i'],single_units[`i']) if rank==`i' & remaining[`i']>0
	quietly drop remaining
	quietly egen sumNew=sum(new_num_systems)
	quietly gen remaining=sumsystems-sumNew
	quietly drop sumNew
}

gen newTotalDamages=new_num_systems*total_damages

egen SumNewtotal_damages=sum(newTotalDamages)

***BOOTSTRAP Damages
foreach var of varlist total_damages1b1-total_damages8b99 {
	gen boot_total_damages_run = `var'*new_num_systems
	egen SumBoot`var' = sum(boot_total_damages_run)
	gen gain_`var' = SumBoot`var' - SumOld`var'
	drop SumBoot`var' SumOld`var'
	drop boot_total_damages_run
}
egen gain_p25 = rowpctile( gain_total_damages1b1- gain_total_damages8b99), p(2.5)
egen gain_p975 = rowpctile( gain_total_damages1b1- gain_total_damages8b99), p(97.5)
egen gain_mean = rowmean( gain_total_damages1b1- gain_total_damages8b99)
egen gain_sd = rowsd(gain_total_damages1b1- gain_total_damages8b99)
gen gain_pt_est = SumNewtotal_damages - sumOldPtEst
drop gain_total_damages1b1- gain_total_damages8b99

putexcel set `"${data_final}/tableA5.xlsx"', sheet(sheet1) modify

putexcel B2 = gain_mean/1000000
putexcel C2 = gain_sd/1000000
putexcel D2 = gain_p25/1000000
putexcel E2 = gain_p975/1000000